Optimal Combinations of Control Strategies and Cost-Effectiveness Analysis of Dynamics of Endemic Malaria Transmission Model

In this study, we propose and analyze a determinastic nonlinear system of ordinary differential equation model for endemic malaria disease transmission and optimal combinations of control strategies with cost effective analysis. Basic properties of the model, existence of disease-free and endemic equilibrium points, and basic reproduction number of the model are derived and analyzed. From this analysis, we conclude that if the basic reproduction number is less than unity, then the disease-free equilibrium point is both locally and globally asymptotically stable. The endemic equilibrium will also exist if the basic reproduction number is greater than unity. Moreover, existence and necessary condition for forward bifurcation is derived and established. Furthermore, optimal combinations of time-dependent control measures are incorporated to the model. By using Pontryagin's maximum principal theory, we derived the necessary conditions of optimal control. Numerical simulations were conducted to confirm our analytical results. Our findings were that malaria disease may be controlled well with strict application of the combination of prevention of drug resistance, insecticide-treated net (ITN), indoor residual spray (IRS), and active treatment. The use of a combination of insecticide-treated net, indoor residual spray, and active treatment is the most optimal cost-effective and efficacious strategy.


Introduction
Malaria is an old deadly parasitic disease and transmitted to the human through the bites of infected female Anopheles mosquito. It is still treated as a serious challenge to the global world population. According to 2020 World Health Organization (WHO) report, there were 241 million cases and 627 thousand deaths from malaria globally. Specifically, according to the report, 80% of the total deaths occurred in Africa among children aged under 5 years old due to malaria [1].
The most popularly used malaria intervention strategies in sub-Saharan Africa and Asia include the use of insecticide-treated bed net (ITN) and indoor residual spray (IRS) against the vector and treatment of infected individuals with drugs. However, despite intensive control efforts, global incidence of malaria is increasing in malaria endemic area. This is due to one of the factors contributed by the spread and emergence of drug-resistant parasite strains against the most common and affordable antimalarial drugs and resistance of malaria vectors against insecticides in sub-Saharan Africa [2,3].
Mathematical modeling has become an important tool in understanding the complex dynamics of disease transmission and in decision-making processes regarding intervention programs for disease control. Concerning malaria disease, Ross developed the first mathematical model. He focused his study on mosquito control and showed that for the disease to be eliminated, the mosquito population should be brought below a certain threshold [4]. Later, the idea of Ross is extended by McDonald to account for super infection [5]. Ngwa and Shu are formulated and analyzed a mathematical model for endemic malaria with variable human and mosquito populations [6]. Wedajo et al. derived and analyzed a deterministic model for the inclusion of infected immigrants on the spread and dynamics of malaria transmission [7]. Chiyaka et al. derived and analyzed the effects of treatment and drug resistance on the transmission dynamics of malaria in endemic areas [8]. Tumwiine et al. analyzed a mathematical model for the transmission and spread of drug-sensitive and drug-resistant malaria strains within human populations [9]. Other studies are carried out by using optimal control theory. Okosun et al. derived and analyzed a malaria disease transmission mathematical model that includes treatment and vaccination with waning immunity and applied optimal control to study the impact of a possible vaccination with treatment strategies in controlling the spread of malaria [10]. Agusto and Khan derived and analyzed optimal control strategies for dengue transmission [11]. Okosun and Makinde modeled the impact of drug resistance in malaria transmission and its optimal control analysis [12]. Bonyah et al. present "Modeling the effects of heavy alcohol consumption on the transmission dynamics of gonorrhea with optimal" [13]. Makinde and Okosun applied optimal control to study the impact of chemotherapy on malaria disease with infective immigrants [14]. Okosun et al. applied optimal control strategies and cost-effectiveness analysis of a malaria model [15]. Keno et al. derived and analyzed optimal control and cost effectiveness analysis of SIRS malaria disease model with temperature variability [16].
In this paper, we study SITRS-SI and SIRS-SI endemic malaria transmission model with standard incidence law that was presented by [12]. Furthermore, we modified the model [6] by omitting the incubating class from the system and incorporating four time-dependent control measures, the class infective in treatment individuals and infectious classes with drug-sensitive and drug-resistant individuals. The purpose of this study is to (i) investigate the stability for both disease-free equilibrium and endemic equilibrium (ii) develop effective ways for controlling the malaria disease (iii) explore the best strategy in terms of reducing the number of malaria infectious populations to zero and minimizing costs

Model Description and Formulation
To study the transmission and spread of malaria in two interacting population of humans (the host) and mosquitoes (the vector), we formulate a model which subdivides the total human population size at time t, denoted by N h ðtÞ into susceptible S h , infected with drug-sensitive malaria strain I hs , infected with drug-resistant malaria strain I hr , infective in treatment T h , and recovered R h . Hence, we have N h t ð Þ = S h t ð Þ + I hs t ð Þ + I hr t ð Þ + T h t ð Þ + R h t ð Þ: Unlike human population, we divide the mosquito population into two compartments: susceptible S v and infected I v .Thus, the total size of the mosquito population at any time t is denoted by Note that S h = S h ðtÞ, I hs = I hs ðtÞ, I hr = I hr ðtÞ, T h = T h ðtÞ, R h = R h ðtÞ, S v = S v ðtÞ, I v = I v ðtÞ, N v = N v ðtÞ, and N h = N h ðtÞ, where ε denotes the rate of infective in treatment individuals that enter in to the susceptible human compartment due to the administered drug kills off the parasites in his/her body and θ represents the rate of recovered individuals that enter in to the susceptible human compartment due to loss of immunity. The susceptible human compartment class is increased due to (i) birth or immigration with a constant requirement rate Λ h ; (ii) the administered drug kills off the parasites in the cell of malaria-infected individuals, and these are found in the infective in treatment human compartment class with respective rate ε; and (iii) immunity loss of individuals in the recovered human compartment class. This group is further reduced by natural death rate μ h and per capita rate of force of infection λ h . Hence, the rate of change of the susceptible human compartment class is given by Furthermore, the rate of change of I hs class is increased by ρð1 − πÞ λ h S h and decreased by the quantity ðμ h + δ h + α + γ s ÞI hs , where μ h , δ h , γ s , and α are natural death rate, diseaseinduced death rate, recovery rate, and the proportion at which I hs enjoy T h , respectively. π is the modification parameter, and ρ is the proportionality per capita rate of susceptible humans entering I hs . Thus, the rate of change is given by Furthermore, the rate of change of I hr class is increased by ð1 − ð1 − πÞρÞλ h S h and decreased by the quantity ðμ h + δ h + γ r + σÞI hr , where μ h , δ h , γ s , and σ are the natural death rate, disease-induced death rate, recovery rate, and proportion at which I hr enjoy T h , respectively. ð1 − ρÞ is the proportionality per capita rate of susceptible humans entering I hr : Thus, the rate of change is given by Also, T h class is increased by αI hs and σI hr but decreased by the quantity ðδ h + μ h + εÞT h , where π 1 is the modification parameter and ε is treatment rate. Thus, the rate change of state equation is given by The rate of change of the population of the partially immune group (recovered class) R h is increased by γ s I hs and γ r I hr and decreased by the quantity ðθ + μ h ÞR h , where 2 Computational and Mathematical Methods in Medicine θ is per capita rate of immunity loss. Thus, the rate change of state equation is given by The rate of change of the population of the susceptible mosquito is generated by the recruitment of mosquitos Λ v and decreased by the natural death rate μ v S v and force of infection λ v S v between the susceptible and infected classes such that Finally, the rate of change of the population of the infected mosquito is generated by the force of infection λ v S v and decreased by the quantity ðμ v + δ v ÞI v , where μ v and δ v are the natural and disease-induced death rates of the mosquito, respectively. Thus, it is given by The notation where β h is the rate of probability of human getting infected, ϕ is the mosquito contact rate with human, and ω is mosquito biting rate.
where β v is the probability of a mosquito getting infected and σ h is the modification parameter. Assumptions in the model build-up are as follows: We represent the above model description diagrammatically in Figure 1. Figure 1 can be expressed or written in seven systems of ordinary differential equations as follows: with initial conditions I hs 0 ð Þ = I 0hs ,
Proof. The total human population at time t is given by N h = S h + I hs + I hr + T h + R h .
After differentiating both sides of N h , 3 Computational and Mathematical Methods in Medicine which gives In the absence of mortality due to malaria, equation (15) becomes Equation (16) is equivalent to The resulting differential inequality can be solved by separation of variables to give Taking the initial conditions t = 0 and denoting N h ð0Þ by N 0h , then the complete solution is given by As t ⟶ ∞, On the other hand, if N 0h > Λ h /μ h , then From (20)- (22), we have Thus, the total human population is bounded in the region: The total mosquito population at time t is given by After differentiating both sides of N v , which gives In the absence of mortality due to malaria, equation (26) becomes Equation (27) is equivalent to The resulting differential inequality can be solved by separation of variables to give Taking the initial conditions t = 0 and denoting N v ð0Þ by N 0v , then the complete solution is given by Computational and Mathematical Methods in Medicine As t ⟶ ∞, On the other hand, if From equations (31)-(33), we have Thus, the total mosquito population is bounded in the region: Thus, the feasible solution set of the system equation of the model enters and remains in the region: Therefore, the basic model (12) is well posed epidemiologically and mathematically. Hence, it is sufficient to study the dynamics of the basic model in Ω.

Positivity of Solutions.
In this section, we aim to obtain nonnegative solutions when dealing with human populations. Therefore, the next discussion below targets on the conditions under which the model being studied has nonnegative solutions.
Proof. To show the positivity of solutions, it is enough to show that each of the trajectories of system (12) is nonnegative for all t ≥ 0: Since the right-hand side of system (12) is completely continuous and locally Lipschitzian on C 1 , the solution ðS h ðtÞ, I hs ðtÞ, I hr ðtÞ, T h ðtÞ, R h ðtÞ, S v ðtÞ, I v ðtÞÞ of system (12) with initial condition equation (13) exists and is unique on [0, k), where 0 < k < +∞.
It follows from equation (3) that the differential inequality describing the evolution of the susceptible human population over time t is given by Hence, Thus, From equation (4), we have ðdI h Þ/dt ≥ −ðμ h + δ h + α + γ s ÞI hs ðtÞ that is equivalent to I hs ðtÞ ≥ exp ½−ðμ h + δ h + α + γ s Þt > 0.
From equation (5), From equation (6), From equation (7), Thus, From equation (9), Therefore, we can see that S h ðtÞ > 0, I hs ðtÞ > 0, I hr ðtÞ > 0, T h ðtÞ > 0, R h ðtÞ > 0, S v ðtÞ > 0, and I v ðtÞ > 0 for all t ≥ 0: 3.3. Disease-Free Equilibrium (DFE). In this section, we establish the existence of the disease-free equilibrium point of system (12). The disease-free equilibrium point of the model is its steady-state solutions without infection or disease. Let the right-hand side of equation (12) equal to zero, that is, Then, system (12) becomes Again letting I hs = I hr = T h = R h = I v = 0, then (44) leads to Then by rearranging equations (45) and (46) and after substituting each other, we got Then, the disease-free equilibrium point is given by 3.4. The Basic Reproductive Number ðR 0 Þ. In this section, we will determine the threshold parameter that governs the spread of the disease which is called basic reproduction number. To obtain the basic reproduction number R 0 of model (12), we used the next-generation matrix method so that it is the spectral radius of the nextgeneration matrix [17,18]. The model equations are rewritten starting with newly infective classes: Let X i = I hs I hr I v ð Þ T , i = 1, 2, 3, and X 1 = I hs , X 2 = I hr , and X 3 = I v ; then, by the principle of the nextgeneration matrix, we can obtain The new infection matrix F is given by Computational and Mathematical Methods in Medicine and the transition matrix V is given by Thus, The eigenvalue of FV −1 can be obtained.
Then, the eigenvalues are From λ 1, , λ 2 , and λ 3 , the dominant eigenvalue is λ 2 . Therefore, the basic reproductive number is given by

Local Stability of Disease-Free Equilibrium
Theorem 3. The disease-free equilibrium point F 0 of system (12) is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1.
Proof. The local stability of the system is determined by the signs of the eigenvalues, and it is further proved by linearizing to obtain its Jacobian at disease-free steadystate points so that the Jacobian matrix of system (12) at disease-free equilibrium point F 0 is defined and given by 7 Computational and Mathematical Methods in Medicine 2   6  6  6  6  6  6  6  6  6  6  6  6  6  6  4   3   7  7  7  7  7  7  7  7  7  7  7  7  7 where To obtain the eigenvalue of (57), we compute Det ðFV −1 − λIÞ = 0 as follows: As the first and sixth columns contain only the diagonal terms which form the two negative eigenvalues, −μ h and −μ v , the other four eigenvalues can be obtained from the following: In the same way, the third and fourth columns contain only the diagonal terms which form the two negative eigenvalues, −k 3 and −k 4 , and the other four eigenvalues can be obtained from the following: Thus, the eigenvalues of the submatrix are the roots of the characteristic equation Here, By the principle of Routh-Hurwitz criteria [19], equation (61) has negative real eigenvalues if and only if b > 0, d > 0, and bc > d. Clearly, we see that b > 0 because it is the sum of positive variables, but d > 0 if and only if 1 − R 2 0 > 0 which is equivalent to R 0 < 1, and hence, all the determinant of eigenvalues of equation (57) will have negative real eigenvalues. Therefore, the disease-free equilibrium point F 0 is locally asymptotically stable.

Global Stability of a Disease-Free Equilibrium (DFE)
Proof. To prove Theorem 4, we consider the Lyapunov function. So let Then by differentiating L, both sides lead to After substituting ðdI hs Þ/dt,ðdI hr Þ/dt, and ðdI v Þ/dt from (12) into (64) and simplifying it, then we get Since Since Computational and Mathematical Methods in Medicine So, dV/dt ≤ 0 if ðR 2 0 − 1Þ ≤ 0 which leads to R 0 ≤ 1. dV/ dt = 0 if and only if I hs = I hr = 0 or R 2 0 = 1: Therefore, by LaSalle's invariant principle [20], every solution to equations of model system (12) with initial condition (13) in Ω approaches to the disease-free equilibrium point F 0 at time t leading to infinity whenever R 0 ≤ 1: and it occurs when the disease persists in the community. To obtain it, we equate all the model equations in (12) to zero. Then, we obtain the following: Substituting (68) into (70), we get Substuting (68) and (72), respectively, into (69), we get where R 0 is the basic reproduction number given by (4.2) and Equation (72) admits a trivial solution λ * h = 0 which corresponds to the disease-free equilibrium point (DFEP). Now, we assume that λ * h ≠ 0, and the existence of endemic equilibria is regulated by the quadratic equation (73) It follows that the number of endemic equilibria of (12) depends on the coefficients b 0 , b 1 , and b 2 as follows.
It follows here also that when we put the value of Λ h = 0:071 and μ h = 0:115 and use Table 1 for the values of other parameters, the two roots are presented graphically as shown in Figure 2.
From Figure 2, we note that there is stable disease-free region when λ * h = 0, and when R 0 = 1, the force of infection mosquitoes to humans starts to increase in stable endemic region where also the disease starts to spread again and hence forward bifurcation. Here, the blue line represents stable equilibrium, and the red line represents unstable equilibrium.

Existence of Bifurcation.
To show the existence of bifurcation of (12), we used the method developed by [21][22][23][24]. In these approaches, the center manifold theory was employed 9 Computational and Mathematical Methods in Medicine to determine the existence and types of bifurcation using the coefficients, a and b, of the normal form representing the dynamics of the system. These coefficients decide the bifurcation. In particular, if a < 0 and b > 0, then the bifurcation is forward; if a > 0 and b > 0, then the bifurcation is backward.
Using this approach, the following result may be obtained. Theorem 6. Model equation (12) exhibits a forward bifurcation at R 0 = 1 if the condition B 0 < 0 is satisfied. Here, B 0 is given by Proof. Let us choose the rate of transmission of infection from an infectious mosquito to a susceptible human with β h as the bifurcation parameter. We observe that R 0 = 1 is equivalent to and the linearized Jacobian matrix evaluated at F 0 and β * h is denoted and given by Here, Þare the solutions of the characteristic equation Det½ J F 0 , β * h ð Þ− λI = 0. Thus, the determinant is given by If and only if λ 1 = −μ h < 0, λ 2 = −μ v < 0, λ 3 = −k 2 < 0, and λ 4 = −k 4 < 0, where c 0 = 1, If we substitute 1 (one) for R 0 into equation (83) or (84), then system (12) will have a simple zero eigenvalue, and the  (12) can then be written in the form dx/dt = FðxÞ; here, F = f 1 , f 2 , f 3 , f 4 , f 5 , f 6 , f 7 À Á T as shown below: The coefficients a and b are defined in Theorem Appendix B.1. of [21]. That is, To compute the coefficient equations a and b, we determine the right and left eigenvectors corresponding to the zero eigenvalue. Thus, the components of the right eigenvectors denoted by w i , for i = 1, ⋯, 7, are given by Equation (87) can be written as Solving equation (88) gives Similarly, the components of the left eigenvector denoted by v i , for i = 1, ⋯, 7, are given by Equation (90) can be written as Solving equation (91) gives Taking into account system (85) and considering only the nonzero components of the left eigenvector v, it follows that the f is denotes the right-hand side of system (85). It can be checked that Now, using the nonzero components of left eigenvector v, right eigenvector w, and the nonzero partial derivatives h ÞðF 0 ÞÞ, the coefficients a and b are determined as Here, B 0 denotes the following parametric expression: Clearly, it can be observed that b is absolutely a positive quantity while a is negative if B 0 < 0. Hence, it is a forward bifurcation.
Therefore, we conclude that the model equations in (12) exhibit forward bifurcation at R 0 = 1 if the condition B 0 < 0 is satisfied.

Computational and Mathematical Methods in Medicine
The Routh-Hurwitz criteria for polynomial equation (97) will give six negative eigenvalues if the conditions given below are satisfied: e i > 0, for i = 1, 2, 3, ⋯, 6. The relevant Routh-Hurwitz criteria in [19,25] could be used to show that model system (12) is locally asymptotically stable when R 0 > 1:

Extension of the Model into an
Optimal Control 4.1. Optimization Process. For malaria, control efforts are carried out to limit the spread of the disease and free the community from the burden of the disease. Optimal control theory is a method that has been widely used to solve for an extremum value of an objective function involving dynamic variables. In this section, optimal control theory of [26] has been applied in deriving optimal control problem that incorporates the timedependent control measures, namely, prevention of drug resistance, insecticide-treated bed net (ITN), treatment, and indoor residual spray (IRS) as functions of time to investigate the optimal control efforts needed to control malaria disease. The controls are assigned reasonable upper and lower bounds and practiced on time interval [t 0 , t f ], where t 0 and t f are the initial and final time, respectively. The control variables as used in the basic malaria model (12) are described as follows: The preventive control measure drug resistance at time t is denoted by u 1 ; this is aimed at minimizing the proportion of the emergence of drug resistance of malaria strains as well as spread of the disease dynamics. It can be implemented by improving the way drugs are used through improving prescribing, follow-up practices, and patient compliance.
The preventive control measure insecticide-treated bed net (ITN) at time t is denoted by u 2 ; it has twofold effects on the malaria disease dynamics: (i) reduces the number of bites from mosquitoes as they physically provide a barrier between the infectious mosquitoes and (ii) reduces the population of the mosquitoes by killing them after they land on the treated net.
The control measure treatment with drugs at time t is denoted by u 3 ; this includes treating individuals who developed symptoms of the disease.
The preventive control measure indoor residual spray (IRS) at time t is denoted by u 4 as preventive measure; i.e., insecticide spray on the breeding site of mosquitoes reduces the number of mosquito populations by killing the rest indoors after feeding.
After incorporating the above controls into basic model system (12), we get the following modified state equation: Here, we also use following the approach developed by [26] and proposed the following objective function J which is used to minimize the number of infected humans with drug-sensitive and drug-resistant malaria parasite strains, infective in treated human populations, and total mosquito populations while keeping the costs of applying the controls u 1 , u 2 , u 3 , and u 4 as low as possible.
where i = 1, 2, 3, 4 and A 1 , A 2 , A 3 , and A 4 and d 1 , d 2 , d 3 , and, d 4 are the coefficients associated to the state variable and controls, respectively. Following the approach of [26], the cost of the controls has been chosen to be quadratic. Thus, the goal is to find an optimal control quadruple, u * 1 , u * 2 , u * 3 , and u * 4 such that where U = f u 1 ðtÞu 2 ðtÞu 3 ðtÞ u 4 ðtÞ : 0 ≤ u i < 1, i = 1, 2, ⋯,

Computational and Mathematical Methods in Medicine
Pontryagin's maximum principle [27] converts system (99) with (100) and (101) into a problem of minimizing pointwise the Hamiltonian H with respect to u 1 , u 2 , u 3 , and u 4 : Thus, the Hamiltonian H equation is given by Here, LðI hs , I hr , T h , I v , u 1 , u 2 , u 3 , u 4 , tÞ = A 1 I hs + A 2 I hr + and λ i , for i = 1, 2, 3, 4, 5, 6, 7, are adjoint variable. Using the existence result for the optimal control [28,29], we established the following theorem as follows.

29
Computational and Mathematical Methods in Medicine numerically by applying Runge-Kutta fourth-order schemes of the approach [26]. The implementation of the scheme was done using the MATLAB package.
The parameter values provided in Table 1  In order to estimate the impact of the incorporated timedependent control measures on the malaria disease dynamics, we proposed four optimal combinations of control strategies. These strategies have been chosen based on their effectiveness in minimizing the spread of the malaria disease. To find the solution, we applied numerical methods, that is, systems (99) and (103) with (105) numerically solved by using the parameter values given in Table 1 and MATLAB software and applying Runge-Kutta fourth-order schemes. Thus, the proposed optimal combinations are as follows: (1) Strategy a: combination of the use of preventive control of drug resistance, insecticide-treated net (ITN), and treatment of infective individuals Strategy a: control with the prevention of drug resistance, insecticide-treated net (ITN), and treatment (u 1 ≠ 0, u 2 ≠ 0, u 3 ≠ 0, and u 4 = 0) In this strategy, we compare the strategy a situation where no control (u 1 = 0, u 2 = 0, u 3 = 0, and u 4 = 0) was used with the application of strategy a. It can be seen from (g) Figure 6: Simulations of the model showing the effect of preventive control of drug resistance, insecticide-treated net, indoor residual spray, and treatment controls.  Figures 3(a)-3(f) that there is a significant increase in the number of susceptible and recovered human populations and a significant decrease in the number of infected with drug-sensitive strains, infected with drug-resistant strains, infective in treated human populations, and infected mosquito populations compared to the strategy with no control at a given time, respectively. From this, one can observe that strict application of strategy a for a period between 10 and 30 days is sufficient to reduce the number of individuals with malaria symptoms and malaria-infected vectors to zero. It can be noted that a combination of the prevention of drug resistance, insecticide-treated nets ITN, and treatment can play an important role in minimizing malaria infectious. The control profile shown in Figure 3(g) shows that controls u 1 , u 2 , and u 3 decrease from the maximum of 100% to the lower bound. This suggests that a high effort is required for preventive control of drug resistance u 1 , insecticidetreated net (ITN) u 2 , and medical treatment u 3 of individuals under this strategy. Strategy b: control with the prevention of drug resistance, indoor residual spray (IRS), and treatment (u 1 ≠ 0, u 2 = 0, u 3 ≠ 0, and u 4 ≠ 0) In this strategy, we compare the strategy a situation where no control (u 1 = 0, u 2 = 0, u 3 = 0, and u 4 = 0) was used with the application of strategy b. It can be seen from Figures 4(a)-4(f) that there is a significant increase in the number of susceptible and recovered human populations and a significant decrease in the number of infected with drug-sensitive strains, infected with drug-resistant strains, infective in treated human populations, and infected mosquito populations compared to the strategy with no control at a given time, respectively. Even though this strategy min-imizes the number of malaria infectious populations, however, it is not enough to eliminate the disease at a given time, and hence, there is a need for additional control effort to eliminate the disease out of the community. The control profile shown in Figure 4(g) shows that controls u 1 , u 3 , and u 4 decrease from the maximum of 100% to the lower bound. This suggests that a high effort is required for preventive control of drug resistance u 1 , indoor residual spray (IRS) u 4 , and medical treatment u 3 of individuals under this strategy.
Strategy c: control with insecticide-treated net (ITN), indoor residual spray (IRS), and treatment (u 1 = 0, u 3 ≠ 0, u 3 ≠ 0, and u 4 ≠ 0) In this strategy, we compare the strategy a situation where no control (u 1 = 0, u 2 = 0, u 3 = 0, and u 4 = 0) was used with the application of strategy c. It can be seen from Figures 5(a)-5(f) that there is a significant increase in the number of susceptible and recovered human populations and a dramatic decrease in the number of infected with drug-sensitive strains, infected with drug-resistant strains, infective in treated human populations, and infected mosquito populations compared to the strategy with no control at a given time, respectively. With the application of strategy d, I v within 8 days, I hs within 10 days, T h within 11 days, and I hr within 30 days will be eliminated from the system. This result is a bit more promising than strategy a and strategy b. The control profile shown in Figure 5(g) shows that control u 3 is at 55% initially and decreases from the maximum of 65% to the lower bound while controls u 2 and u 4 decrease from the maximum of 100% to the lower bound within 90 days. This suggests that a high effort is required for the use of insecticide-treated net u 2 and indoor residual   Strategy d: control with the prevention of drug resistance, insecticide-treated net (ITN), indoor residual spray (IRS), and treatment (u 1 ≠ 0, u 2 ≠ 0, u 3 ≠ 0, and u 4 ≠ 0) In this strategy, we compare the strategy a situation where no control (u 1 = 0, u 2 = 0, u 3 = 0, and u 4 = 0) was used with the application of strategy d. It can be seen from Figures 6(a)-6(f) that there is a significant increase in the number of susceptible and recovered human populations and a significant decrease in the number of infected with drug-sensitive strains, infected with drug-resistant strains, infective in treated human populations, and infected mosquito populations compared to the strategy with no control at a given time, respectively. With the application of strategy d, I v , I hs , T h and I hr within timet = 8,10,11 and 30days respectively will be eliminated from the system. This result is a bit more promising than when strategy a and strategy b except possibly strategy c which yielded almost the same results. The control profile shown in Figure 6(g) shows that controls u 1 , u 2 , u 3 , and u 4 decrease from the maximum of 100% to the lower bound. This suggest that a high effort is required for preventive control of drug resistance u 1 , insecticide-treated net u 2 , indoor residual spray (IRS) u 4 , and medical treatment u 3 of individuals under this strategy.

Cost-Effectiveness Analysis.
To analyze the costeffectiveness of the strategies, we employ the approach of incremental cost-effective ratio (ICER) in [15,41]. The ICER is used to compare the cost and the health outcomes of two alternative intervention strategies that compete for the same resources.
Based on the values of ICER of the strategy k, the alternative that is more expensive and less ineffective is then excluded. This is done after simulating the optimal control model and then ranking strategies in order of increasing effectiveness measured as the total infection averted. We calculate ICER based on strategy k ðk = a, b, c d Þ by using the following formula: The total infection averted by optimal strategy k during the time period t f is denoted and defined as where And A kI hs , A kI hs , A kI hs , and A kI v are the total infected humans with drug-sensitive disease parasite strains, infected humans with drug-resistant disease parasite strains, infective in treated human populations, and infected mosquito populations, respectively, averted by optimal strategy k during the time period final time t f , I * * hs , I * * hr , T * * h , and I * * v that are the optimal solutions associated with optimal controls ðu * 1 , u * 2 , u * 3 , u * 4 Þ (49), and I hs ð0Þ, I hr ð0Þ, T h ð0Þ, and I v ð0Þ are the corresponding initial values. Note that these initial values are obtained as the equilibrium proportion of system (49) with no postexposure intervention ðu 1 = u 2 = u 3 = u 4 = 0Þ. The total cost associated with a strategy is denoted and given by where d 1 corresponds to the per person unit cost following preventive control of drug resistance intervention, d 2 corresponds to the per person unit cost of preventive control of ITN intervention, d 3 corresponds to the per person unit cost of treatment intervention, and d 4 corresponds to the per person unit cost of preventive control of IRS intervention. Parameter values from Table 1 are used to estimate the total infection averted and total cost presented in Table 2.
It follows here that after rearranging control strategies in Table 2 in increasing order of effectiveness ðTA k Þ, incremental ð113Þ Table 3 incorporates all the four strategies, that is, strategies a, b, c, and d, and corresponding total infection averted during each strategy or effectiveness ðTA k Þ from Table 2 in increasing order and incremental cost-effectiveness ratio ICER.
Because strategy a is less effective and more costly among the other strategy, we excluded it in Table 4. Thus, Table 4 incorporates only three strategies (that is, strategies b, c, and d) and corresponding total infection averted during each strategy or effectiveness ðTA k Þ from Table 2 in increasing order and incremental cost effectiveness ratio ICER.
Similarly, because strategy d is less effective and more costly among strategies b, c, and d, we excluded it from Table 5. Thus, Table 5 incorporates only two strategies, that is, strategies b and c, and corresponding total infection averted during each strategy or effectiveness ðTA k Þ from Table 2 in increasing order and incremental cost effectiveness ratio ICER.
Finally, we conclude that the comparison result reveals that strategy c is cheaper than strategy b. Therefore, strategy c (the use of preventive control of insecticide-treated net (ITN), treatment for infected with drug-sensitive parasite, infective in treatment individuals, and indoor residual spray (IRS)) is the best of all possible strategies due to its more effectiveness and less cost or due to its cost-effectiveness and healthy benefits.

Discussions and Conclusions
In this study, a nonlinear system of ordinary differential equation model that describes the dynamics of malaria disease transmission is formulated and analyzed. Basic properties of the model, existence of disease-free and endemic equilibrium points, and basic reproduction number of the model are derived and analyzed. From this analysis, we conclude that the basic reproduction number of the model is the threshold parameter between the extinction and persistence of the disease. That is, if the basic reproduction number is less than unity, then the disease-free equilibrium point is both locally and globally asymptotically stable resulting in the disease removing out of the host populations. The endemic equilibrium will also exist if the basic reproduction number is greater than unity. Moreover, existence and necessary condition for forward bifurcation when the basic reproduction number is equal to unity are derived and established. Further-more, optimal combinations of time-dependent control strategies, that is, prevention of drug resistance, insecticide-treated nets ITN, treatment, and indoor residual spray (IRS), are incorporated to the model. Using optimal control theory of Pontryagins's maximum principle, existence, necessary conditions, and optimality system for the optimal control problem are derived and analyzed.
The optimality system, that is, state system (99) and adjoint system (103) with optimal controls in (105), was solved numerically by applying Runge-Kutta fourth-order schemes. Numerical simulation results of these show that the use of a combination of prevention of drug resistance, insecticide-treated net (ITN), active treatment, and indoor residual spray (IRS) or strategy d performs the best by controlling the disease for the time period of intervention introduced. Moreover, cost-effectiveness analysis using ICER indicates that strategy c is the most optimal cost-effective and efficacious strategy.
Finally, we conclude that the application of optimal combinations of control strategies not only reduces the number populations with malaria symptoms but also reduces the emergence of drug resistant malaria strains as well as the spread of the disease.

Data Availability
In this paper, some parameter values are cited from the previous literatures which are listed in our references. Because of the lack of the available literatures and data, the other parameter values are technically assumed (see details in Table 1).